#using Pkg
#Pkg.status()
using LinearAlgebra
using Plots
using DiffEqDevTools,BenchmarkTools
using OrdinaryDiffEq
using JLD2, FileIO
using Parameters,NBInclude
include("../../IRKGaussLegendre.jl/src/IRKGaussLegendre.jl")
using .IRKGaussLegendre ## bertsio lokala exekutatzeko
PATH_SRC="../src/"
include(string(PATH_SRC,"MyBenchmarksTools.jl"))
include(string(PATH_SRC,"IRKGL_SIMD.jl"))
using .IRKGL_SIMD
WARNING: using IRKGL_SIMD.launch_method_tests in module Main conflicts with an existing identifier. WARNING: using IRKGL_SIMD.WPTests in module Main conflicts with an existing identifier. WARNING: using IRKGL_SIMD.launch_IRKGL_simd_tests in module Main conflicts with an existing identifier.
function RigidBody!(du,u,p,t)
I1=p[1]
I2=p[2]
I3=p[3]
du[1] = I1*u[2]*u[3]
du[2] = I2*u[1]*u[3]
du[3] = I3*u[1]*u[2] + 0.25*sin(t)^2
end
t0=0.
tF=100.0 # 10.0
tspan=(t0,tF)
tspan_B=(BigFloat(t0),BigFloat(tF))
p = [-2.0,1.25,-0.5]
u0 = [1.0;0.0;0.9]
p_B = BigFloat.([-2.0,1.25,-0.5])
u0_B = BigFloat.([1.0;0.0;0.9])
dim=length(size(u0))
udim=u0[1:1]
prob = ODEProblem(RigidBody!,u0,tspan,p)
prob_B= ODEProblem(RigidBody!,u0_B,tspan_B,p_B)
sol = solve(prob,Vern9(),abstol=1/10^14,reltol=1/10^14);
plot(sol, xlabel="t", title="Rigid Body problem")
#sol =solve(prob_B,Vern9(),save_everystep=false, abstol=1e-30,reltol=1e-30)
#@save "./Data/Rigid_Body_sol100.jld2" sol
#test_sol = TestSolution(sol)
@load "./Data/Rigid_Body_sol100.jld2" sol
test_sol = TestSolution(sol)
final_state=sol.u[end]
(Float32(sol.t[end]))
100.0f0
tols=abstols=reltols=1.0 ./ 10.0 .^ (6:14)
dts8= 4*[0.01, 0.015, 0.02, 0.025, 0.03, 0.04, 0.05]
dtsVern=dts8/3;
nruns=1000
s=8
wp1=launch_IRKGL16_tests(IRKGL16(),final_state, prob, dts8; nruns=nruns)
wp2=launch_IRKGL_simd_tests(final_state, prob, dim, s, dts8,
initial_interp=-1, partitioned=false, nruns=nruns)
wp11=launch_method_tests(Vern9(), final_state, prob, tols, nruns=nruns)
wp12=launch_method_tests(Vern9(), final_state, prob, dtsVern, adaptive=false, nruns=nruns);
yrange=(10^(-3.0),10^(-1.75))
plot(title="Rigid Body ",xlabel="Error", ylabel="Time (s)", ylims=yrange,
xticks=[1e-16,1e-15,1e-14,1e-13,1e-12, 1e-11, 1e-10])
#
plot!(wp1.errors, wp1.times, seriestype=:scatter, scale=:log10, label="", color="blue")
plot!(wp1.errors, wp1.times, scale=:log10, label="", lw=3, color="blue") #label="IRKGL16_seq"
#
plot!(wp2.errors, wp2.times, seriestype=:scatter, scale=:log10, label="",color="green")
plot!(wp2.errors, wp2.times, scale=:log10, label="", lw=3, color="green") #label="IRKGL16_simd"
#
#
plot(title="Rigid Body",xlabel="Error", ylabel="Time (s)")
#
plot!(wp11.errors, wp11.times, seriestype=:scatter, scale=:log10, label="Vern9-Adap",color="red")
plot!(wp11.errors, wp11.times, scale=:log10, label="", lw=3, color="red")
#
plot!(wp12.errors, wp12.times, seriestype=:scatter, scale=:log10, label="Vern9-Fix",color="blue")
plot!(wp12.errors, wp12.times, scale=:log10, label="", lw=3, color="blue")
#
#
plot!(wp2.errors, wp2.times, seriestype=:scatter, scale=:log10, label="IRKGL_simd_s8)",color="green")
plot!(wp2.errors, wp2.times, scale=:log10, label="", lw=3, color="green")
#
yrange=(10^(-3.5),10^(-2.25))
plot(title="Rigid Body",xlabel="Error", ylabel="Time (s)", ylims=yrange,
xticks=[1e-16,1e-15,1e-14,1e-13,1e-12, 1e-11, 1e-10])
#
#plot!(wp11.errors, wp11.times, seriestype=:scatter, scale=:log10, label="Vern9-Adap",color="red")
#plot!(wp11.errors, wp11.times, scale=:log10, label="", lw=3, color="red")
#
plot!(wp12.errors, wp12.times, seriestype=:scatter, scale=:log10, label="",color="orange") #label="Vern9"
plot!(wp12.errors, wp12.times, scale=:log10, label="", lw=3, color="orange")
#
#
plot!(wp2.errors, wp2.times, seriestype=:scatter, scale=:log10, label="",color="green") #label="IRKGL16_simd",
plot!(wp2.errors, wp2.times, scale=:log10, label="", lw=3, color="green")
#
#